Diluted 3d-Random Field Ising Model at zero temperature with metastable dynamics 
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We study the influence of vacancy concentration on the behaviour of the three dimensional Ran- 
dom Field Ising model with metastable dynamics. We focus our analysis on the number of spanning 
avalanches which allows for a clean determination of the critical line where the hysteresis loops 
change from continuous to discontinuous. By a detailed finite size scaling analysis we determine the 
phase diagram and estimate numerically the critical exponents along the whole critical line. Finally 
we discuss the origin of the curvature of the critical line at high vacancy concentration. 
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I. INTRODUCTION 

Externally driven systems at low enough temperature 
often display rate independent hysteresis. This out-of- 
cquilibrium phenomenon occurs because intrinsic disor- 
der creates multiple energy barriers that the system can- 
not overcome due to the very weak thermal fluctuations. 

The study of zero temperature models with metastable 
dynamics has been very succesful for understanding rate 
independent hysteresis. A prototype case is the Random 
Field Ising Model (RFIM) with single spin-flip relaxation 
dynamics pi Q- Although the model is formulated in 
terms of magnetic variables (external field H and mag- 
netization m) it can be applied to the study of many phe- 
nomena associated to low temperature first-order phase 
transitions in disordered systems, e.g. martensitic trans- 
formations JU, fluid adsorption in porous solids Q, fer- 
roelectrics |5(, etc. 

Disorder is a intriguing concept: in the RFIM it is 
introduced via independent and quenched random fields 
on each lattice site, gaussian distributed with zero mean 
and standard deviation a. In real materials, disorder 
is much more complicated and includes features at all 
length scales: vacancies, interstitials, composition fluctu- 
ations, dislocations, strain fields, grain boundaries, sam- 
ple surfaces, edges and corners, etc. Thus it is interesting 
to add to the RFIM other sources of disorder in order to 
see how the non-equilibrium behaviour is modified. 

The goal of this paper is to study the diluted RFIM 
at T — with metastable dynamics and analyze the con- 
sequence of introducing a concentration c of quenched 
vacancies. The interplay between the two kinds of disor- 
der (random fields and vacancies) will be at the origin of 
the properties of the a — c phase diagram. 

One of the striking results concerning the RFIM with 
metastable dynamics, as already pointed out in the sem- 
inal paper of Sethna et al.Q, is the occurence of a crit- 
ical point when the amount of disorder a is increased. 
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The to vs. H hysteresis loops change from discontinuous 
(like in a ferromagnet) when a < a c to continuous (like 
in a spin-glass) when a > a c . This result was demon- 
strated using mean-field analysis and numerical simula- 
tions in 3d systems. This problem was also studied within 
the Renormalization Group formalism 0, QJ • Moreover, 
many properties of the critical point have been also stud- 
ied analytically on Bethe lattices |llllTfl[llllT^|. Exper- 
imental evidence for the occurence of such a critical point 
has been found in different magnetic systems j 1 ■' 3. 1 1 4| . 

Another interesting result of the RFIM with 
metastable dynamics is that it reproduces the experimen- 
tal observation that the m(H) trajectories of such ather- 
mal systems are discontinuous at small scales. The evo- 
lution proceeds by avalanches from a mestastablc state to 
another. In the RFIM the avalanche size distribution be- 
comes a power- law at the critical point. Experimentally, 
scale free distributions of avalanche properties have been 
found in many systems HE EE [13, HM O E3, E3| 
A first attempt to study the influence of dilution in such 
avalanche size distributions was done some years ago [2j| . 
The results of this work, however, should be considered 
as only qualitative, given the fact that the studied sys- 
tem was two-dimensional |28| . the analysis focused only 
on the avalanche distributions and the results concerning 
the phase diagram were very approximate. 

The order parameter that vanishes at the critical point 
is the size of the macroscopic discontinuity Am. The 
analysis of this quantity from simulations is very intrin- 
cate. In finite-size systems it is very difficult to make 
the distinction between a macroscopic jump and a mi- 
croscopic avalanche. The measured order parameter only 
displays reasonable finite-size scaling (FSS) properties 
when the simulated systems are very large [24). Recent 
studies [2^, H(| , have shown how the critical point can be 
characterized in systems of moderate size. The key point 
is to detect the so-called "spanning" avalanches which are 
the magnetization jumps that involve a set of spins that 
spans the whole finite system (e.g. cubic lattice) from 
one face to the opposite one. By this method avalanches 
in finite systems can be classified as non-spanning, ID- 
spanning, 2D-spanning, or 3D-spanning. The average 
numbers N\, N2 of ID and 2D-spanning avalanches dis- 
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play a peak at a value of a that shifts with system size L 
and tends to a c when L — ► oo. The numerical data can 
then be scaled according to the FSS hypothesis pEj 



N a = L e N a (uL 1/v ) 



(1) 



where a = 1, 2. The exponent v = 1.2 ±0.1 characterizes 
the divergence of the correlation length (£ ~ (cr — cr c ) _i/ ) 
whereas = 0.10 ± 0.02 characterizes the divergence of 
the number of critical avalanches. The scaling variable 
u(a) is analytic and measures the distance to the critical 
point. It can be fitted by the second order expression 



u(a) = 



A 



(2) 



with cr c = 2.21 and A = -0.2. The behaviour of the 3D- 
spanning avalanches is more complex because they are of 
two different kinds: (i) critical 3D-spanning avalanches 
that behave like the ID and 2D ones and (ii) subcritical 
3D-spanning avalanches that will corespond to the Am 
discontinuity in the thermodynamic limit. The analysis 
is more difficult and requires a double finite-size-scaling 
technique. This will not be used in the present paper. 
Instead we will focus only on the behaviour of the average 
numbers iVi and JV2 in the presence of vacancies and 
propose a FSS hypothesis by using a bivariate scaling 
variable u(a, c) that allows to study the full a— c diagram. 

In section ITT1 we define the model and the dynamics. 
In section IIIII we present results of the numerical simu- 
lations. In section llVl we formulate the FSS hypothesis 
and determine the critical line. In section we propose 
approximations to the bivariate scaling variable u(a, c). 
In section IVil we discuss the interplay between vacancies 
and avalanches and finally, in section TVIII we summarize 
our main findings and conclude. 



II. MODEL AND SIMULATIONS 

The diluted 3D-RFIM on a cubic lattice with N sites 
(JV = L x L x L) is defined by the following Hamiltonian 
(magnetic enthalpy): 



H 



ciCjSiSj 



JV 



N 



hiCiSi - -ff Sji 



(3) 



where Si = ±1 are Ising spin variables, Cj = 0,1 indi- 
cates the presence of a vacancy (cj = 0) or not (c, = 1) 
on each site, hi are quenched random fields gaussian dis- 
tributed with zero mean and standard deviation a and H 
is the driving field. The first sum extends over all distinct 
nearest-neighbour (n.n.) pairs. Vacancies are quenched 
and randomly distributed over the lattice. Their concen- 
tration is measured by c = 1 — J^i C i/N. 

The metastable dynamics is implemented as follows: 
the system is externally driven by the field H which is 
adiabatically swept from — 00 where the system is fully 



negatively magnetized (Si = —1) to +00. (Si = +1). 
The spins flip according to a local relaxation dynamical 
rule, 
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H) 



(4) 



where the sum extends over all the n.n. of Si- When a 
spin flips, it may trigger an avalanche. The unstable spins 
are flipped synchronously until a new stable situation is 
reached. 

The hysteresis loop is obtained by computing the mag- 
netization 



iV 



= Y^S lCl /N 



(5) 



as a function of the applied field H. Magnetization 
avalanches are recorded along the whole increasing field 
branch and their spanning properties are analyzed by 
using "mask" vectors (as explained in Ref. |25|) that al- 
low to classify them as non-spanning, ID-spanning, 2D- 
spanning and 3D-spanning. In this work we shall mainly 
study the number of spanning avalanches of each kind 
which are recorded in the full upwards branch. These 
numbers, N±, N% and N3 which depend on L, a and c 
corresponds to averages over more than 10 realizations 
with different random fields and random vacancy posi- 
tions. The disorder averages are denoted by the symbol 
(•}. We study sytems of sizes ranging from L = 8 to 
L — 64 in a number of points on the a — c diagram, as 
indicated schematically in Fig. ^ 




Region II 



FIG. 1: (Color online) Coordinates of the points studied by 
numerical simulations in the a — c diagram. The finite size 
scaling analysis presented in section llVl is performed in regions 
I and II. 



III. NUMERICAL RESULTS 

The general evolution of the average hysteresis loops as 
a function of a and c is shown in Fig. El Oue can observe 
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the transition from discontinuous loops to smooth loops 
when tr or c are increased. It can also be seen that the 
saturation magnetization decreases with increasing c. 
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FIG. 2: (Color online) Average hysteresis loops correspond- 
ing to a system of size L = 32 for different values of a and c 
as indicated. 



Figure |3 shows the behaviour of the coercive held (H coe ) 
as a function of the concentration of vacancies for differ- 
ent values of a. As can be seen, (H coe ) decreases with 
increasing c and increasing a. The behaviour with in- 
creasing c exhibits an inflection point at the transition as 
can be seen in the inset of Fig. [3] which shows the numer- 
ical derivative of the (H coe ) with respect to c. Such an 
inflection point does not exist in the non-diluted model 
when the coercive field is plotted as a function of a. This 
feature, which can be of interest for the determination 
of the critical point in experiments, is probably related 
to the fact that (H coe ), as a function of c should van- 
ish at c < 1 whereas as a function of a it only vanishes 
assymptotically when a — > oo. 



Fig. Q] shows the distribution D(s;cr,c,L) of avalanche 
sizes (the size s of an avalanche is the number of spins 
flipped) for the same cases as in Fig. |2 in log-log scale. 
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FIG. 3: (Color online) Coercive field as a function of the 
vacancy concentration c for different values of the amount of 
disorder a. The inset shows the behaviour of the numerical 
derivative dH coe /dc which exhibits a maximum on the tran- 
sition line. Data correspond to averages in a system of size 
L — 64. 

The histograms include all avalanches irrespective of 
their spanning properties. The qualitative picture is that 
power law distributions are obtained along a critical line 
with an exponent that seems to be the same for all val- 
ues of c. Apparently, no differences can be observed when 
comparing the transition induced by changing a from the 
transition induced by changing c. Below the critical line 
the distributions show a peak for large values of s which 
correspond to the ID, 2D and 3D spanning avalanches. 
Above, the distributions have an exponentially damped 
character. 



Figure shows the average number of ID, 2D and 3D 
spanning avalanches as a function of a for increasing val- 
ues of the vacancy concentration c ranging from to 0.5. 
Data corresponds to a system with size L = 16. 



The same information is displayed in Fig. [S]for a system 
with size L = 48. 



The behaviour for small and intermediate vacancy con- 
centration is qualitively similar to that found for the non- 
diluted model [25j . The average numbers Ni and N2 dis- 
play peaks, whereas N$ shows a peak on the edge of a 
step function. Note that for L = 16 the peak height in 
Ni(a, c, L) and N 2 (cr, c, L) seem to decrease with increas- 
ing c. This behaviour, however, is much less apparent for 
larger systems (L = 48). Thus, it is possibly due to a fi- 
nite size effect. 

At higher concentrations (c > 0.4) Aq and JV 2 be- 
gin to develop a flat plateau at low a. The reason for 
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FIG. 4: (Color online) Avalanche size distributions corre- 
sponding to a system with L — 32 at different values of c and 
a as indicated. Data are represented in log-log scale. 



this plateau can be well understood by looking at the 
3d-plot in Fig. which represents the average number 
Ni(a,c,L) for L — 32. The plateau in the constant c 
cuts of Figs. [S]and[|)]is due to the fact that the crest of 
the iVi and N2 functions bends towards the a = axis. 



IV. FINITE SIZE SCALING HYPOTHESIS 

The hypothesis that we want to check numerically is 
that, in the presence of vacancies, the critical point found 
at c = transforms into a critical line for a wide range of 
concentrations. Thus the critical exponents found previ- 
ously should be equally valid for the description of the 
behaviour of the average numbers N\ and N2 with c > 0. 
According to this hypothesis we shall propose the follow- 
ing corresponding FSS behaviour: 



N a (a,c,L) = L e N a {uL 1/v ) 



(0) 



where a = 1, 2 and u(cr, c) is a bivariate scaling variable 
measuring the distance to the critical line. The expo- 
nents 6 and v as well as the functions N a were already 




FIG. 5: (Color online) Average number of ID-spanning 
avalanches (a), 2D-spanning avalanches (b) and 3D-spanning 
avalanches (c) as a function of a for different values of the va- 
cancy concentration c, as indicated by the legend. Lines are 
guides to the eye. Data corresponds to numerical simulations 
of a system with size L = 16. 



found in previous works 25] . Therefore, the hypothesis 
is quite strong and indicates that all the A^i and N2 data 
corresponding to different sizes L, different vacancy con- 
centrations c and different amounts of disorder a must 
collapse into an already known function. The only free- 
dom is in the determination of the bivariate scaling vari- 
able u that should be analytic. Before constructing it 
in the next section, we can make a first test of Eq. El by 
cheking the scaling on the critical line. Note that, setting 
u = 0, Eq. becomes: 



N a (a,c,L)=L N a (O) 



(7) 



where a and c should be on the critical line. Since we 
know (from Ref. |H that 9 = 0.10, iVi(0) = 0.12 and 
7V 2 (0) = 0.07, we can deduce that the different curves 
N a (a 7 c, L)/L e N a (0) should cross at height 1 on the crit- 
ical line, independently of L. Two examples are shown 
in Fig. [H] that corresponds to two cuts (one at constant 
a and the other at constant c) in the a — c diagram. As 
can be seen, the critical line can be determined with high 
accuracy. By analyzing a large number of such a and c 




FIG. 6: (Color online) Average number of ID-spanning 
avalanches (a), 2D-spanning avalanches (b) and 3D-spanning 
avalanches (c) as a function of a for different values of the va- 
cancy concentration c, as indicated by the legend. Lines are 
guides to the eye. Data corresponds to numerical simulations 
of a system with size L = 48. 



FIG. 7: (Color online) Surface plot representing Ni(a,c,L) 
for L — 32. The dashed lines on the basal plane represent the 
position of the cuts in Fig|H]at c = 0.15 and a = 0.9. 



For small values of a, the critical line displays a vertical 
behaviour. The critical value of the vacancy concentra- 
tion c c above which the hysteresis loops do not display a 
discontinuity can be fitted to c c = 0.426 ± 0.003 



V. BIVARIATE SCALING VARIABLE 

In general the bivariate scaling variable is a function 
that can be expanded as: 

u(a, c) = ao + a\o + a 2 c + a^ac + a^a 2 + a^c 2 + . . . (8) 



cuts we have constructed it. The result is shown in Fig. 

Note that the process can be repeated independently 
with N\ and N 2 . The two independent lines overlap al- 
most perfectly. 



The obtained critical line is linear up to c ~ 0.3. A 
least squares fit gives <r c (c) = <7 C (0) + Ac with <r c (0) = 
2.21 ±0.01 and A = -4.09 ±0.03. The value er c (0) = 2.21 
is in total agreement with the previous estimate for the 
non-diluted model [25|. 

It is remarkable that the finite size scaling hypothesis 
allows the collapse of the data up to large values of c, 
far from the point c — where the scaling function and 
the exponents where determined. It is also remarkable 
that scaling works even after the bend that is observed 
for c > 0.3. (Note that the crossing point shown in Fig. 
EI a) corresponds to a value of a where the critical line is 
not linear.) 



Since c and a are not necessarily very small along the 
critical line, it is difficult to know a priori how many 
terms in the expansion will be needed in order to find a 
good scaling collapse. The direct determination of a large 
number of coefficients from the numerical data is difficult. 
Therefore we shall take a different strategy taking into 
account, as much as possible, the previously known data. 

As a first step we will concentrate in the region c < 0.3 
where the coexistence line shows a linear behaviour and 
we will try to use an expansion up to quadratic terms 
only. By forcing that the condition u = is satisfied on 
the fitted coexistence line, we deduce that u satisfies: 

u((T, c) = (a — (T c — Xc)(b + b\o + b 2 c) (9) 

We should also consider the fact that the scaling variable 
is known to be well described by a second order expansion 
(up to a 2 ) for c = as indicated in Eq. [21 After some 
algebra one can determine the two parameters bo and b\ : 

b = (l-A)/a c = 0.543 ±0.002 (10) 
bi = A/a 2 = -0.041 ± 0.001 (11) 
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FIG. 8: (Color online) Examples of crossing points on the 
critical line on cuts (a) parallel to the c axis and (b) parallel 
to the a axis. The different symbols correspond to different 
system sizes as indicated by the legend. Continuous lines are 
guides to the eye. The horizontal dashed-line indicates the 
height 1 where the curves are supposed to cross according to 
Eq.0 
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FIG. 9: (Color online) Critical line in the a — c diagram de- 
termined from the crossing points in Ni (•) and 7^2 (x). The 
dashed line is the fit discussed in the text, and the thin dis- 
continuous lines indicate the cuts along which the correlation 
in Fig. E)is computed. 
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Therefore we are left with a single free parameter hi that 
should allow to collapse into a single curve all the data in 
the scaling region, for different values of er, c and L. We 
have considered all the available data in region I of Fig. 
n Note also that the same hi parameter must be used to 
scale both Ni and iV2 data. The best two collapses are 
shown in Fig. EH and Fig. ITTlfor b 2 = -0.13. 



FIG. 10: (Color online) Finite-size-scaling collapse of the 
average number of ID-spanning avalanches in region I. The 
continuous line shows the Lorentzian function in Eg. 1121 



on top of the data in Figs. HOIandllll 
Ni(x) ' 



(1.73 - 0.53a; + 0.10a; 2 ) 3 ' 



(121 



Note that the data for c = are also included in this plot. 
Therefore we are obtaining two scaling functions N\ and 
N~2 compatible with those of Ref. |23- In that reference 
the scaling functions were approximated by Gaussians, 
although it was also shown that there were systematic 
deviations. In this work we have tried to fit the data with 
more complex functions (with three free parameters). We 
have found a very good x 2 by using the following modified 
Lorentzians, which are represented by a continuous line 



1 



(1.83 - 0.59a; + 0.13a; 2 ) 



2^4.6 



(13) 



As a second step we try to build u(a, c) for the data very 
close to the a = axis. In this region II (see Fig. |5J the 
transition line is again quite linear and in fact is almost 
vertical. This means that to measure the distance to 
the critical line it should be sufficient to use the variable 
(c — c c ). We have considered the following second-order 
expansion: 



u(c) 



D 



(14) 
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FIG. 11: (Color online) Finite-size-scaling collapse of the 
average number of 2D-spanning avalanches in region I. The 
continuous line shows the Lorentzian function in Eg. 1131 



FIG. 13: (Color online) Finite size scaling collapse of the 
average number of 2D-spanning avalanches in region II. The 
continuous line shows the lorentzian funcion in Eg. 1131 



Note that k' is not a free parameter. It can be fixed 
by imposing that the definitions of the scaling variables 
© and {TIJ coincide at a — and c = 0. Thus k' = 
(A— 1)/(B — 1). The only free parameter for the collapse 
of the data is B. Best collapses are shown in Fig. IT21 
and 1 131 for N\ and N2 respectivelly using the best choice 
B = -0.2 (thus k' = 1.) 
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FIG. 12: (Color online) Finite size scaling collapse of the 
average number of ID-spanning avalanches in region II. The 
continuous line shows the lorentzian funcion in Eq. 1121 



The continuous lines in both figures correspond to the 
same lines as in Figs. 1101 and 1111 We can thus conclude 
that we have built two good approximations (given by 
Eq. IrjlandfHI) in regions I and II, to the unique scaling 
variable which will display a more complex behaviour 
in the intermediate region were the critical line bends, 
probably with higher order terms. 



VI. DISCUSSION 

In this section we try to understand why the critical 
line exhibits such a curvature. There must be a physical 
reason that goes beyond the mere effect of the dilution of 
the system and unstabilizes even more the phase with 
the ferromagnetic-like discontinuity. We propose that 
the effect is related to the percolation of vacancies above 
c. p = 0.3116, a value which is, indeed, very close to the 
limit where the critical line loses its linearity. To justify 
this hypothesis numerically we have studied, for each par- 
ticular realization of disorder the distribution of the clus- 
ters of vacancies and the position of the avalanches. In 
particular we have determined the spatial position of the 
largest vacancy cluster (which above c p will correspond 
to the percolating cluster in the thermodynamic limit). 
It is clear that the neigbouring sites of this percolating 
cluster of vacancies are an easy path for the propagation 
of an avalanche since these sites have a smaller number 
of neighbours. To distinghuish such sites we have defined 
a local flag that takes values 6^ = 1 when a site belongs 
to the border of the largest cluster of vacancies or 6j = 
otherwise. We have also recorded the largest avalanche 
during the i7-scan (which will correspond to the spanning 
avalanche below the critical line in the thermodynamic 
limit) and marked its position with a flag ej = 1. With 
these two variables we have defined the correlation be- 
tween the border of the largest cluster of vacancies and 
the largest avalanche as: 

(££eifr«)-(££ei)(££iM 

sfii £ t) ~ ii £ e *> 2 v £ ^ - £ b ^ 2 

(15) 



Note that since £j and bi take values 1, only, the power 
2 in the first bracket inside the square roots can be sup- 
pressed. This correlation is equal to 1 when the spanning 
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FIG. 14: (Color online) Correlation between the border of 
the largest cluster of vacancies and the largest avalanche as a 
function of c. Data corresponding to different system sizes are 
represented by different symbols as indicated by the legend. 
The curves correspond to cuts in the phase diagram at a — 1 
(a), a = 0.5 (b) and a = 0.1 (c). 



avalanche sits exactly on the border of the spanning clus- 
ter of vacancies. The behaviour of p^^ as a function of 
c is shown in Fig. ^] for three different values of a that 
correspond to the dashed lines indicated in Fig. Inland for 
increasing system sizes as indicated by the legend. The 
important observation is that the curves for a — 0.5 and 
(7 = 0.1 exhibit two crossing points. One is located at c p 
and the other on the critical line (it thus shifts with a). 
For a concentration of vacancies below c p or above the 
critical line, the behaviour of the curves with increasing 
L indicates that the correlation vanishes in the thermo- 
dynamic limit, whereas in the region between the two 
crossing points the correlation increases with increasing 



system size. A value p = 1 is probably not reached since 
the spanning avalanche is larger than the border of the 
percolating cluster of vacancies. By this analysis, we have 
thus identified the origin of the curvature of the critical 
line: when vacancies percolate the spanning avalanche 
propagates along the border of the percolating cluster 
of vacancies. The propagation in such a constrained en- 
vironment decreases the amount of disorder needed to 
break the infinite macroscopic avalanche into small mi- 
croscopic jumps. However, as shown in the previous sec- 
tion this mechanism does not changes the values of the 
critical exponents. 



VII. SUMMARY AND CONCLUSIONS 

We have analyzed the influence of dilution on the criti- 
cal properties of the 3D-RFIM at T = with metastable 
dynamics. We have shown that the critical point associ- 
ated to the change in the shape of the hysteresis loop from 
discontinuous to continuous loops becomes a critical line 
which we have located in the a — c phase diagram. The 
critical properties close to this line are characterized by 
the same critical exponents as in the non-diluted model. 
This result indicates that it should be possible to find RG 
arguments showing that there is a unique fixed point at 
T = in the disorder parameter space.that includes, at 
least, both random fields and dilution |2|. We have com- 
puted quadratic approximations to the scaling variable 
in two different zones of the phase diagram that allow 
for a bivariate finite-size-scaling collapse on a universal 
scaling function. Finally we have proposed an explana- 
tion for the the curvature observed in the critical line 
when the concentration of vacancies increases above the 
percolation limit: the spanning avalanche that is respon- 
sible for the discontinuity of the hysteresis loops, has a 
tendency to follow the neigbourhood of the percolating 
cluster of vacancies. 
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